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ABSTRACT 


The  objective  of  this  program  was  to  investigate  the  risk-quantified  design  methods  for  the 
design  and  certification  of  the  structure  for  future  military  aircraft.  The  current  Department  of 
Defense  (DoD)  process  for  developing  and  supporting  aircraft  structures  relies  heavily  on 
testing,  including  full-scale  testing,  to  ensure  structural  integrity  and  reliability.  This  process 
adds  cost  and  time  to  any  new  development  program  or  service  life  extension  effort.  In  the 
future  DoD  vision  of  lean,  agile  acquisition  programs,  a  new  paradigm  is  needed  for  making 
decisions  about  the  airworthiness  of  an  aircraft.  Other  industries,  such  as  nuclear  power  and 
offshore  oil,  have  developed  a  risk  quantification  process  for  deciding  if  a  given  structure  should 
be  placed  in  service.  A  similar  risk  quantification  process  could  be  applied  in  the  airworthiness 
decision  process  for  small  aircraft  fleets  and  hypersonic  aircraft  after  adapting  it  to  the  DoD 
procurement  process,  the  aircraft  design  process,  and  the  safety  issues  particular  to  aircraft. 

The  Risk  Quantified  Structural  Design  and  Evaluation  -  Warm  Primary  Structure  (RQSD&E) 
program  was  conducted  in  two  phases.  The  first  phase  was  exploratory  determining  the 
advantages  of  risk  quantified  structural  evaluations  of  simple  structural  configurations.  The 
second  phase  focused  on  developing  simple  techniques  to  accomplish  the  following  three  things: 

1)  Identify  structural  items  that  affect  vehicle  reliability  or  safety  (important  structural 
items), 

2)  Design  structural  components  for  a  specified  probability  of  failure,  and 

3)  Validate  structural  models  with  a  limited  number  of  data  containing  variability,  or 
uncertainty. 

This  project  determined  that  the  challenge  for  risk-quantified  design  and  evaluation  does  not 
seem  to  be  in  developing  new  more  sophisticated  probabilistic  techniques,  but  in  getting  the 
existing  techniques  into  common  practice.  It  needs  to  become  relatively  easy  for  the  average 
engineer  without  a  lot  of  expertise  in  probabilistic  technology  to  use  probabilistic  methods.  The 
way  to  do  this  is  to  make  incremental  changes  to  existing  design  and  maintenance  processes  by 
enabling  engineers  to  readily  obtain  more  and  better  information  upon  which  to  base  any 
decisions  about  the  airworthiness  of  an  aircraft  upon. 
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1.0  INTRODUCTION 


The  objective  of  this  program  was  to  investigate  the  risk-quantified  design  methods  for  the 
design  and  certification  of  the  structure  for  future  military  aircraft.  The  current  Department  of 
Defense  (DoD)  process  for  developing  and  supporting  aircraft  structures  relies  heavily  on 
testing,  including  full-scale  testing,  to  ensure  structural  integrity  and  reliability.  This  process 
adds  cost  and  time  to  any  new  development  program  or  service  life  extension  effort.  This  has 
been  acceptable  because  the  cost  of  half  a  dozen  test  aircraft  and  the  associated  testing  cost  and 
time  is  insignificant  compared  to  the  several  hundred  aircraft  being  bought.  Compared  to  the 
production  schedule  for  these  large  numbers  of  aircrafts  the  several  years  of  testing  time  was 
acceptable.  However,  the  future  vision  for  DoD  aircraft  purchases  is  lean,  agile  programs  that 
efficiently  produce  less  than  ten  aircraft  of  a  given  type.  In  this  situation,  the  addition  of  a  single 
airframe  for  testing  can  increase  the  program  cost  by  10%  or  more.  And  the  years  required  to 
complete  the  tests  could  be  longer  than  the  production  schedule.  Therefore,  another  process  is 
needed  to  ensure  aircraft  structural  integrity  in  these  lean,  agile  acquisition  programs. 

Other  industries,  such  as  nuclear  power  and  offshore  oil,  offer  alternatives  to  the  current  build- 
and-test  paradigm.  In  both  industries,  structural  integrity  is  an  important  attribute.  Because  of 
the  uniqueness  of  each  structure  or  its  operating  conditions,  testing  is  not  a  viable  approach  for 
evaluating  structural  integrity.  Both  industries  have  developed  a  risk  quantification  process  for 
deciding  if  a  given  structure  should  be  placed  in  service.  A  similar  risk  quantification  process 
could  be  applied  in  the  airworthiness  decision  process  for  small  aircraft  fleets  and  hypersonic 
aircraft  after  adapting  it  to  the  DoD  procurement  process,  the  aircraft  design  process,  and  the 
safety  issues  particular  to  aircraft.  The  risk  quantification  process  must  work  within  the  Systems 
Engineering  Process. 
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2.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 


The  Risk  Quantified  Structural  Design  and  Evaluation  -  Warm  Primary  Structure  (RQSD&E) 
program  was  conducted  in  two  phases.  The  first  phase  was  exploratory  determining  the 
advantages  of  risk  quantified  structural  evaluations  of  simple  structural  configurations.  The 
second  phase  focused  on  developing  simple  techniques  to  accomplish  the  following  three  things: 

1)  Identify  structural  items  that  affect  vehicle  reliability  or  safety  (important  structural 
items), 

2)  Design  structural  components  for  a  specified  probability  of  failure,  and 

3)  Validate  structural  models  with  a  limited  number  of  data  containing  variability,  or 
uncertainty. 

2.1  Exploration 

The  objective  of  this  phase  was  to  demonstrate  the  benefits  of  probabilistic  design  and  risk 
assessment  for  aircraft  structural  design.  Two  benchmark  problems  were  selected  to  demonstrate 
the  benefits:  the  thermal  buckling  response  of  a  rectangular  panel  [1],  and  the  thermal  buckling 
response  of  the  wing  skin  of  a  conceptual  hypersonic  aircraft  [2].  In  addition  to  using  these 
benchmark  problems  to  investigate  the  benefits  of  probabilistic  design  and  risk  assessments, 
practical  methods  for  reliability-based  design  of  aerospace  structures  were  investigated. 

The  exploratory  investigations  and  their  detailed  results  were  reported  in  two  reports  [1,2]. 

These  investigations  concluded  that  probabilistic  analysis  provides  a  way  to  determine  how 
variability  in  loading,  geometry,  materials,  and  environment  affect  the  reliability  of  a  design  and 
the  contribution  of  each  design  parameter.  Probabilistic  design  requires  a  number  of  analyses, 
typically  finite  element  analyses  for  structures.  General  purpose  probabilistic  analysis  software 
packages,  such  as  NESSUS®  [3]  and  UNIPASS™  [4],  that  can  automatically  control  the 
analysis  process  are  commercially  available.  These  codes  also  provide  information  on  the 
sensitivity  of  the  results  to  various  input  parameters.  This  sensitivity  information  is  useful  in 
driving  towards  more  robust  designs. 

There  are  several  major  issues  which  deter  the  use  of  probabilistic  design  practices.  These  issues 
include: 


•  defining  reliability  design  criteria  for  the  system  and  the  individual  components; 

•  handling  geometric  variability  that  requires  new  finite  element  models  and  meshes;  and 

•  determining  the  probability  distributions  for  design  parameters  with  limited  information. 

2.2  Technique  Development 

The  techniques  developed  during  the  second  phase  of  this  program  have  been  reported  in  several 
journal  papers  and  conference  presentations.  Only  those  techniques  not  presented  elsewhere  are 
presented  in  detail.  References  are  provided  to  papers  and  presentations  discussing  the 
remaining  techniques. 
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2.2.1  Identifying  Important  Structural  Items 

Important  Structural  Items  are  defined  as  struetural  details  that  either  eontribute  greatly  to  the 
reliability  of  the  airframe,  Reliability  Critical  Structural  Items  (RCSIs),  or  significantly  impact 
the  structural  integrity  and  safety  of  the  airframe,  Significant  Structural  Items  (SSIs).  As 
depleted  in  Figure  1,  these  two  sets  of  struetural  details  eomprise  two  subsets  of  the  entire 
airframe.  RCSIs  are  the  loeations  and  eomponents  with  a  higher  probability  of  failure  and 
significantly  affect  the  overall  system  reliability.  SSIs  may  or  may  not  have  a  high  probability  of 
failure;  however,  if  an  SSI  fails  the  aircraft  will  be  lost.  Hence,  some  details  may  be  contained  in 
both  subsets,  but  it  is  not  neeessary  for  the  two  subsets  to  be  eoineident. 


Figure  I:  Depiction  of  Relation  between  Reliability  Critical  Structural  Items  and 

Significant  Structural  Items 

The  goal  of  this  effort  was  to  develop  well-defined  processes  for  identifying  these  two  types  of 
structures  using  the  structural  models  to  supplement  experience  and  engineering  judgment  as  to 
what  struetural  eomponents  are  important.  Struetural  eomponents  identified  as  important  ean 
then  be  treated  with  speeial  attention  during  design,  fabrieation,  assembly  and  maintenanee.  For 
instanee,  important  structural  components  can  be  subjected  to  special  analysis  checks  during 
design,  special  inspections  during  fabrication  and  assembly,  and  identified  as  locations  for  health 
monitoring  sensors  in  serviee. 

2.2.I.I  Identifying  Reliability  Critical  Structural  Items 

A  technique  for  rapidly  screening  structural  failure  modes  and  locations  to  determine  those  that 
contribute  most  to  the  reliability  of  the  airframe  was  developed.  The  method  uses  a  formal  error 
metric  to  quantify  the  error  incurred  with  respect  to  the  system  reliability  if  a  particular  location 
is  filtered,  i.e.,  removed  from  the  set  of  locations  whose  probability  of  failures  are  considered  to 
contribute  to  the  system  probability  of  failure.  Locations  with  errors  below  a  user-defined 
probability  of  failure  threshold  (typically  in  the  range  of  10'  to  10'  )  are  filtered;  these  locations 
contribute  little  to  the  overall  system  reliability.  The  method  integrates  naturally  with  a  first- 
order  reliability  method  (FORM)  representation  of  the  failure  modes  and  is  very  fast  (one  million 
limit  states  can  be  filtered  in  2  hours  on  a  desktop  machine).  A  hierarchical  multi-level  approach 
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has  been  developed  to  verify  the  important  limit  states.  The  procedure  for  performing  this 
screening  has  been  demonstrated  for  an  aircraft  wing  [5-6],  an  exhaust  manifold  for  an  internal 
combustion  engine  [6],  and  has  been  verified  using  Monte  Carlo  sampling.  Further  validation 
efforts  will  include  comparison  against  tests  of  a  T-38  longeron. 

2.2. 1.2  Identifying  Significant  Structural  Items 

A  method  for  identifying  SSIs  using  vulnerability  analysis  was  investigated.  A  structure,  or  any 
part  of  it,  is  vulnerable  when  relatively  small  damage  to  the  structure  leads  to  disproportionately 
large  consequences,  such  as  inability  to  carry  limit  load  or  the  structure  becoming  a  mechanism 
[7].  Preliminary  results  in  Penmetsa,  et  al.  [8]  indicated  potential  to  identify  structural  elements, 
or  regions,  that  are  vulnerable  even  when  they  are  not  highly  loaded  or  have  a  low  probability  of 
failure. 

The  challenge  is  that  numerous  analyses  of  the  structure  are  required  considering  failure  at  each 
structural  location  in  turn.  A  technique  for  performing  these  analyses  without  having  to 
recompile  the  global  stiffness  matrix  in  a  finite  element  model  was  developed  to  make  these 
analyses  less  cumbersome.  The  technique  utilizes  the  fact  that  a  failure  in  a  single  location  does 
not  dramatically  change  the  global  stiffness  matrix. 

2.2.2  Designing  to  a  Specified  Probability  of  Failure 

Most  aerospace  engineers  intuitively  understand  that  using  a  larger  factor  of  safety  (uncertainty) 
and/or  margin  of  safety  will  result  in  a  more  reliable  part.  However,  the  actual  reliability  of  the 
part  has  never  quantified  to  determine  how  much  larger  or  lower  these  factors  need  to  be. 

Simple  relationships  between  the  factor  of  safety  (uncertainty),  or  the  margin  of  safety,  and  the 
probability  of  failure  from  static  overstress  were  developed  for  general  distributions  of  load  and 
strength  in  [9-10].  Similarly,  a  method  for  computing  the  probability  of  fracture  for  a  cracked 
structure  subjected  to  a  defined  loading  distribution  using  a  three  chart  nomograph  was 
developed  [11-13]. 

While  developing  these  methods,  it  became  apparent  that  there  was  a  need  to  find  the  distribution 
for  the  product  or  ratio  of  random  variables.  There  were  no  general  methods  for  making  these 
calculations,  so  a  procedure  was  developed. 

2.2.2.I  Finding  Products  and  Quotients  of  Distributions 

Ravi  Penmetsa 
Wright  State  University 

Determining  the  distribution  of  the  product,  or  ratio,  of  random  variables  is  a  common  operation 
in  structural  reliability  analyses.  The  strength  distribution  and  the  cross  section  area  distribution 
must  be  multiplied  together  to  determine  the  distribution  for  the  failure  load  of  a  structural 
element  so  that  the  probability  that  the  applied  loads  exceed  the  failure  load  can  be  computed. 
Or  the  applied  load  distribution  needs  to  be  divided  by  the  cross  sectional  area  distribution  to 
obtain  the  distribution  of  applied  stress. 
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Most  efforts  towards  developing  methods  for  determining  the  distribution  for  the  product,  or 
ratio,  of  random  variables  have  concentrated  on  methods  that  produce  closed-form,  analytical 
Probability  Distribution  Functions  (PDF)  for  products  of  variables  with  particular  types  of 
distributions  [14-21].  Recently,  methods  capable  of  handling  mixed  families  of  distributions  like 
a  product  of  Gamma  and  Weibull  distributions  have  been  developed  [22].  While  these  methods 
are  limited  to  obtaining  a  PDF  for  the  product  of  two  random  variables,  Podolski  [23]  has 
presented  a  technique  to  determine  the  PDF  of  a  product  of  “n”  variables.  However,  this  method 
is  restricted  to  Gamma  distributions. 

There  are  no  methods  for  situations  where  the  variables  are  from  arbitrary  families  of 
distributions.  Many  times  the  distributions  for  the  random  variables  in  structural  reliability 
analysis  are  different  from  those  distributions  for  which  the  closed  form  solutions  are  available. 
Moreover,  it  is  impractical  to  implement  a  general  structural  reliability  analysis  routine  where 
different  solution  approaches  are  required  for  different  factor  distributions. 

The  method  described  here  is  applicable  to  the  product  or  ratio  of  any  number  of  random 
variables  and  has  no  restriction  on  the  type  of  PDFs.  One  of  the  drawbacks  of  this  methodology 
is  that  it  is  restricted  to  only  non-negative  random  variables.  However,  in  structural  applications 
random  variables  representing  geometric  parameters  and  material  strength  variations  are  always 
positive.  Several  examples  are  presented  to  demonstrate  the  proposed  method. 

2.2.2.1.1  Product  and  Ratio  of  Random  Variables 

Let  Y=XiX2  be  the  product  of  two  random  variables  whose  PDFs  are  defined  by/;  and/2.  The 
PDF  of  Yify)  is  obtained  as  follows: 

1 .  Take  the  logarithm  of  both  sides  of  Y=X]X2: 

log  [T]  =  log  [V;V2]  =  log  [V;]  -t  log  [V2]  (1) 

2.  Recast  the  equation  in  terms  of  new  variables  Z,  A;,  and  A2: 

Z  =  log  [T]  =  log  [V;]  -t  log  [V2]  =Ai+A2  (2) 

3.  Determine  the  PDFs  of  A;  and  A2  using  the  chain  rule  for  PDF  transformation.  This  will 
be  discussed  in  the  following  section. 

4.  Determine  the  PDF  of  the  sum,  Z  =  A;  -1-  A2,  using  the  convolution  integral.  This  is  a 
straightforward  approach  for  standard  distributions  like  normal,  Weibull,  etc.  Since  A; 
and  A2  may  not  have  standard  PDFs,  the  PDF  of  the  sum  is  evaluated  using  Fast  Fourier 
Transformation  (FFT)  based  convolution.  Details  of  the  FFT  convolution  process  are 
presented  in  the  section  titled  “Convolution  Using  Fast  Fourier  Transform”. 

5.  Transform  the  PDF  of  Z  using  the  chain  rule  to  determine  the  PDF  of  Y  =  exp[Z]. 

These  steps  are  applicable  to  both  the  product  and  the  ratio  of  any  number  of  independent 
random  variables  with  any  random  distributions.  The  final  PDF  obtained  using  this  approach  will 
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be  a  discrete  set  of  data  points  representing  the  PDF  values  for  the  range  of  Y.  While  this 
approach  does  not  directly  produce  closed  form  equations  for  the  final  PDF,  a  general  method  for 
fitting  a  weighted  sum  of  beta  distributions  can  be  used  to  determine  a  closed  form  equation  for 
the  final  distribution  if  needed.  Many  structural  reliability  analysis  tools  handle  discrete  PDF 
definitions  without  the  need  for  actual  closed  form  equations  [24] . 

2.2.2.1.2  Transformation  of  Random  Variables 


The  expression  for  the  PDF,yA(^),  of  a  variable  A=g(x)  can  be  determined  by  differentiating  its 
Cumulative  Distribution  Function  (CDF),  Fa(A).  [25] 


fA(A)=Af  (A)  = 

dA 


This  can  be  further  simplified  to 

Using  this  chain  rule,  the  PDF  for  A  =  log(x)  is  found  to  be 


A  =  log(x) 


dA 


(3) 


(4) 


(5) 


Equations  (5)  derive  the  PDF  of  log(x)  given  the  PDF  of  the  variable  x.  The  inverse  operation  of 
determining  the  PDF  for  A  =  exp(x)  is 


A  =  exp  (x) 


dx  _  I  _  1 

dA~  A~ 


/a (^) "  ^ ^ (log(^)) 


(6) 


The  last  transformation  will  be  done  on  a  pointwise  representation  of  the  PDF,/ji:(x). 

2.2.2.1.3  Convolution  Using  Fast  Fourier  Transformation 

Fast  Fourier  Transform  (FFT)  based  convolution  is  a  well-developed  technique  and  has  been 
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used  for  many  decades  in  signal  processing.  The  FFT  algorithm  derives  its  efficiency  by 
transforming  the  physical  domain  into  the  frequency  domain.  Convolution  in  the  physical 
domain  is  converted  to  an  algebraic  product  of  two  signals  (PDFs  in  this  case)  in  the  frequency 
domain.  The  FFT-based  convolution  algorithm  for  a  generic  problem  of  the  form  X-aY  is 
outlined  below. 

Step  1:  FFT-based  convolution  is  applicable  for  the  sum  of  random  variables.  Therefore,  X-aY 
needs  to  be  put  in  the  form  X-i-Z.  This  can  be  done  by  introducing  a  variable  Z  =  -  aY.  The  PDF 
of  Z  can  be  determined  using  chain  rule  as  follows 


fz 


-1 

a 


fr 


\  a  J 


(7) 


Step  2:  Once  the  function  is  expressed  as  a  sum  of  two  random  variables  X  and  Z,  their  PDF’s 
need  to  be  discretized  using  a  common  discretization.  This  discretization  enables 
implementation  of  an  efficient  Discrete  FFT  algorithm.  The  discretization  of  the  PDF’s  is 
determined  based  on  the  number  of  points  and  the  bounds  of  the  variables  used  for  the 
convolution.  Most  FFT  algorithms  are  optimized  to  handle  points  as  powers  of  2.  Therefore, 
(=1024)  to  2  (=4096)  points  will  be  used  in  this  paper.  A  convergence  study  can  be  performed 
to  determine  the  number  of  points  required  to  give  the  desired  accuracy.  Once  the  number  of 
points  is  defined,  the  smallest  range  (UpperLimit  -  LowerLimit)  of  X  and  Z  is  selected  to 
determine  the  discretization  increment.  The  smallest  range  is  used  so  that  both  the  distributions 
will  be  at  least  the  selected  number  of  points. 

This  discretization  step  yields  two  vectors  with  different  numbers  of  elements  depending  on  the 
range  of  X  and  Z.  For  example,  consider  a  case  where  A  is  a  random  variable  with  range  [0,  4] 
representing  the  entire  area  under  its  PDF,  and  Z  has  a  range  [3,  9].  If  the  required  number  of 
points  was  determined  to  be  2'°,  then  the  discretization  increment  would  be 

4-0 

A  =  — ^  =  0.0039. 

2"  (8) 
When  the  PDF  of  Z  is  discretized  using  the  same  increment,  the  size  of  the  vector  for  Z  would  be 

5,  =^71^  =  1539  =  2'“ +515;^  2'“  + 2^  (9) 

^  A 

However,  in  order  to  apply  discrete  FFT  algorithm  the  sizes  of  the  vectors  need  to  be  equal. 
Therefore,  7?  zero-valued  elements  are  added  to  the  upper  end  of  the  A-vector  extending  its  range 
from  [0,4]  to  [0,6] . 

With  this  operation  both  the  vectors  have  equal  numbers  of  elements  but  the  number  of  elements 
in  the  vectors  is  not  a  power  of  2.  Two  final  operations  are  performed  that  add  additional  zeros 
until  the  size  of  each  vector  is  equal  to  a  power  of  2,  and  then  double  the  size  of  the  vectors.  The 
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first  operation,  making  the  size  a  power  of  2,  ensures  better  efficiency  since  the  discrete  FFT 
algorithms  are  optimized  to  deal  with  powers  of  2.  The  second  operation,  doubling  the  size  of 
the  vector,  eliminates  the  circular  convolution  issue  of  discrete  FFT.  Both  operations  increase 
the  range  of  the  vectors. 

For  the  above  example,  the  vectors  X[0,  6]  and  Z[3,  9]  are  of  size  2^'^+2^.  Since  the  next  power 
of  2  is  2' \  2^  zero-valued  elements  are  added  to  both  vectors.  The  new  ranges  of  the  two  vectors 
areX[0,  8],andZ[3,  11]. 

Step  3:  Once  the  discrete  vectors  have  been  constructed,  the  next  step  is  to  obtain  the  FFT  of 
each  of  these  vectors.  FFT  routines  are  standard  functions  in  many  mathematical  software 
packages  and  spreadsheets. 

Step  4:  The  FFT  converts  the  convolution  of  two  PDF’s,  X  and  Z,  into  the  product  of  two 
vectors,  A  and  B.  Both  these  vectors  are  multiplied  element  by  element  using  complex  number 
multiplication  for  the  entire  length  of  the  vector  A  to  obtain  the  resultant  vector,  C. 

Step  5:  The  inverse  FFT  of  the  vector  C  will  yield  the  PDF  of  the  initial  function  represented  by 
X-aY. 


Step  6:  The  range  of  the  final  PDF,  X-aY,  is  determined  from  the  ranges  of  X  and  Z.  Since  the 
final  distribution  is  the  sum  of  the  two  variables,  X  and  Z,  the  range  of  X-aY  is  the  interval  found 
by  adding  the  two  ranges.  Based  on  the  values  selected  earlier,  for  Z[0,  8]  and  Z[3,  11],  then  the 
range  for  X-aY  is  [0-1-3,  8-1-II]  =  [3,  19].  It  should  be  noted  that  the  size  of  the  final  vector  is 
equal  to  the  size  of  the  input  vectors  after  doubling  their  length. 

FFT-based  convolution  is  highly  efficient  once  implemented  and  can  be  applied  to  problems 
with  any  PDF  distributions  without  the  need  for  any  approximations  of  the  PDF.  Therefore,  FFT- 
based  convolution  is  ideal  for  implementing  as  a  general  convolution  routine. 

2.2.2.1.4  Numerical  Examples 

The  following  examples  are  selected  to  demonstrate  the  FFT-based  joint  PDF  modeling  process 
for  the  product  and  ratio  of  random  variables.  The  first  example  is  selected  to  demonstrate  the 
accuracy  of  the  presented  numerical  technique  compared  to  analytical  solutions.  The  second 
example  demonstrates  product  of  four  random  variables.  And  the  third  example  demonstrates  the 
combination  of  product  and  ratio  of  variables. 

Example  1: 

In  this  example,  the  PDF  of  a  random  variable,  Z=XY,  is  determined  using  the  FFT-based 
convolution  technique.  Glen,  et  ah,  [26]  have  developed  an  analytical  solution  for  the  situation 
where  X  ~  Uniform[l,2]  and  Y~Uniform[3,4].  These  same  uniform  distributions  are  considered 
here  to  compare  the  numerical  solution  to  the  analytical  solution.  The  analytical  solution  is: 
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(11) 


PDF  of  Z  = 


fz  ~ 


Log[z]- Log[3], 

■  ^og[4]-Log[3], 
3Log[2]  -  Log[z], 


3<z<4 
4  <  z  <  6 
6<  z  <8 


Following  the  steps  outlined  in  the  previous  section,  the  PDF  of  Z  is  determined  from  the  FFT- 
based  convolution  of  log[X]  and  log[Y].  For  this  example,  points  were  used  to  discretize  the 
PDFs  of  the  two  random  variables;  A  is  equal  to  0.00098.  The  intermediate  distributions  used  to 
find  the  PDF  of  Z  are: 


A  =  log(x),  fi  =  log(T) 

fM)  =  e^  for  AG[0,log(2)] 

fs  ^  for  B  G  [log(3),  log;4)  ] 

Figure  2  shows  the  agreement  between  the  numerical  solution  and  the  analytical  solution. 


Figure  2:  Comparison  of  Analytical  and  Numerical  Solutions  for  PDF  of  Z=XY 
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Example  2: 

In  this  example,  product  of  four  random  variables  is  considered  Z=XiX2X3X4.  The  PDF  of  log[Z] 
is  merely  the  convolution  of  the  PDFs  of  log[X7],  log[X2],  log[Xj],  and  log[X4]: 

log  Z  =  log  X]  +  log  X2  +  log  X3  +  log  X4  (12) 

The  distributions  were  selected  as  follows  X/~Normal[100,10],  X2~Weibull[100,  3.5], 
X?~Rayleigh[25],  Z4~Uniform[75,  100]  to  demonstrate  the  applicability  of  the  method  to  any 
combination  of  distributions.  Since  an  analytical  solution  in  not  available  for  this  situation, 
Monte  Carlo  simulation  was  used  to  estimate  the  PDF  of  the  product. 


The  distributions  used  in  the  convolution  are: 

A  =  logZj,  B  =  \ogX2,  C  =  logZ3,  D  =  \ogX^, 


fAA)- 


1 


loV^ 


exp 


A-1 

2 


^e^-100^" 

10 


f,(B) 

fc(c) 


3.5 

— 


100 

1 

25- 


f  pB  A 

3.5 

3.55- 

1  C 

tiOOy 

-exp 


2C  — 
2 


A  ^c\ 


v25y 


fo{D)  = 


e 


D 


25  ■ 


(13) 


It  can  be  seen  in  Figure  3  that  FFT-based  approach  for  determining  the  PDF  of  the  product  of 
four  arbitrary  distributions  matches  the  PDF  resulting  from  the  Monte  Carlo  simulations. 


10 


Figure  3:  PDF  of  Product  of  Four  Variables 

Example  3: 

In  this  example,  a  combination  of  both  product  and  ratios  of  random  variables  is  considered: 


3 

XjXf  ’  (14) 

log  Z  =  log  Xj  +  log  ^2  -  log  X3  -  2.5  log  X4 . 

Monte  Carlo  simulation  is  once  again  used  to  estimate  the  PDF  of  Z  in  order  to  verify  the  FFT- 
based  approach.  The  following  distributions  were  used  in  this  example:  X7~Weibull[1000,  3.5], 
X2~Rayleigh[25],  X5~Normal[500,  80],  X4~Uniform[25,  50].  The  intermediate  distributions 
needed  for  the  convolution  are: 
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A  =  logXj,  5  =  logZ2,  C=logX3,  D  =  logZ 


4  ’ 


/a  (^)  =  77^:^  exp 


1000 


3.5 

3.5A- 

^1000 ; 

^exp 


25-- 

2 


f 


v25y 


c-i 

2 


^ -100^" 
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(15) 


x,xf 

Figure  4:  PDF  of  a  Combination  of  Product  and  Ratio  of  Random  Variables 


Figure  4  shows  that  the  two  results  match  accurately  over  the  entire  range  of  Z. 
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2.2.2.1.5  Summary  of  Finding  Products  and  Quotients  of  Distributions 

An  FFT-based  convolution  method  for  determining  the  PDF  of  the  product  and  ratio  of  random 
variables  was  found  to  be  applicable  to  any  distribution  of  random  variables.  This  approach  is 
applicable  only  to  non-negative  random  variables.  In  most  structural  applications,  the  random 
variables  are  positive  structural  dimensions  or  material  properties  thereby  making  this  approach 
useful  for  structural  reliability  analysis.  This  approach  can  be  used  along  with  regression  models 
that  have  interaction  terms  that  are  products  of  random  variables.  Therefore,  there  are  numerous 
applications  for  which  this  method  can  be  used. 


2.2.3  Quantifying  Uncertainties  in  Large-Scale  Numerical  Simulations  as  Input  to  a 
Validation  Framework 

Jason  Gruenwald,  Daniel  Nevill  and  Mark  Brandyberry 
University  of  Illinois,  Urbana-Champaign 


2.2.3.I  Introduction 

A  major  obstacle  to  the  use  of  large  scale  multi-physics  simulations  in  design  is  proving  that  the 
simulations  are  reliable  and  valid.  A  methodology  is  needed  to  account  for  the  uncertain  inputs  to 
the  computational  models,  which  cause  difficulty  in  comparing  results  between  the 
computational  models  and  any  available  experimental  data  on  the  physical  system  being 
modeled.  In  order  to  capture  the  variability,  probability  distributions  are  used  so  that  statistical 
measures  of  confidence  can  be  used  on  the  models  to  ensure  they  are  reliable.  Then,  these 
simulations  can  be  used  to  predict  the  responses  from  extreme  conditions  that  would  be  nearly 
impossible  to  reproduce  in  a  lab  setting.  However,  it  is  often  not  feasible  to  directly  produce 
uncertainty  estimates  for  large  computational  models,  since  the  simulations  are  potentially  long- 
running  and  computationally  intensive.  This  project  was  focused  on  developing  a  method  for 
estimating  this  uncertainty  with  a  minimum  of  computational  effort. 

The  goal  was  to  develop  a  large  scale,  multi-physics  methodology  that  quantifies  these 
uncertainties  while  minimizing  the  amount  of  computer  simulation  time.  The  end  result  is  a 
probability  distribution  describing  the  likelihood  of  a  set  of  response  variables  having  certain 
values  (described  by  probability  distributions),  which  could  play  an  integral  part  in  design 
testing.  This  project  investigated  the  use  of  a  clustering  methodology  to  produce  the  uncertainty 
quantification.  The  simulated  system  contained  uncertainty  in  both  fluid  and  structural  elements, 
and  thus  the  Computational  Fluid  Dynamics  (CFD)  modules  from  a  simulation  suite  known  as 
Rocstar  was  coupled  with  the  finite  element  solver  Abaqus  to  model  the  various  aspects  of  the 
problem. 
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Figure  5:  The  internal  structure  of  the  wing 
2.23.2  Description  of  Prohlem 

The  wing  of  a  supersonic  trainer  aircraft  was  chosen  for  demonstration  of  the  process.  The  wing 
has  a  sweep  of  31.9°  and  a  taper  ratio  of  0.2.  It  has  a  constant  NACA  65A004.8  airfoil  from  the 
root  to  tip.  For  the  simulations,  only  the  aerodynamic  and  structural  effects  of  the  wing  at  0° 
flaps  were  considered.  The  aerodynamic  effects  of  the  fuselage  were  ignored.  The  internal 
structure  is  shown  in  Figure  4. 

2.2.3.3  Multi-Physics  Models 

2.2.3.3.1  Fluid  Dynamics 

To  capture  the  uncertainty  in  the  fluid  portion  of  the  simulation,  a  detailed,  three-dimensional 
CFD  package  was  used  to  model  the  air  flow  around  the  wing.  The  Rocstar  suite  is  fully 
integrated  to  enable  three-dimensional,  multi -physics  simulations.  Rocstar  is  comprised  of  two 
CFD  solvers,  Rocflo  and  Rocflu.  Rocflo  is  an  explicit  block-structured  hexahedral  solver,  and 
Rocflu  is  an  implicit  unstructured  mixed-mesh  solver.  For  our  purposes  of  determining  the  lift 
generated  on  the  wing  only  the  structured  block  fluids  solver,  Rocflo,  was  used. 

In  order  to  run  a  CFD  analysis  of  the  airflow  about  the  wing  of  the  aircraft,  a  mesh  was  created 
around  the  area  of  interest.  This  mesh  is  made  of  elements,  and  the  governing  equations  of  flow 
will  be  applied  at  each  element.  The  mesh  created  for  our  purposes  is  fairly  coarse,  in  order  to 
limit  the  total  computation  time  for  this  test  problem.  A  finer  mesh  would  produce  a  more 
accurate  simulation  of  the  flow,  but  would  also  increase  the  computation  time.  The  mesh  used  in 
this  simulation  was  comprised  of  approximately  600,000  elements.  Due  to  symmetry  only  one 
half  of  the  airplane  was  modeled,  and  it  was  assumed  that  the  other  half  is  symmetric  in 
response.  The  grid  was  non-uniform,  body-fitted,  and  the  elements  were  smallest  near  the  wing 
and  increased  in  size  further  away.  This  gave  the  most  accurate  results  at  the  elements  nearest 
the  wing,  where  we  were  most  interested  in  the  flow  characteristics. 
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2.2.3.3.2  Solid  Mechanics 


A  full  scale  idealized  model  of  the  wing,  shown  in  Figure  6,  was  provided  from  the  United  States 
Air  Force.  The  model  was  originally  created  in  NASTRAN,  but  was  translated  into  Abaqus  for 
our  simulations.  It  consisted  of  the  wing  and  the  main  spars  that  attach  the  wing  to  the  fuselage. 
It  modeled  the  structurally  relevant  components,  but  did  not  model  the  leading  edge  or  the 
control  surfaces  on  the  wing.  The  internal  structure  was  idealized  by  beam,  truss,  shell  and 
membrane  elements.  The  wing  was  modeled  with  cantilever  boundary  conditions  with  extra 
pinned  boundary  conditions  along  the  wing  root  to  simulate  the  wing  attachment  to  the  fuselage. 
The  loads  were  generated  through  CFD  simulation  as  introduced  in  section  2. 2. 3. 3.1,  and  were 
applied  as  pressure  loads  through  an  external  transfer  algorithm. 


Figure  6:  The  idealized  model  of  the  wing  was  used  in  Abaqus  for  the  full  scale  simulations 
2.2.3.4  Uncertain  Variables 

The  first  step  in  the  methodology  was  to  determine  which  parameters  were  uncertain  and  how 
they  varied.  Note  that  for  purposes  of  this  project,  the  values  and  distributions  used  needed  to  be 
representative,  but  not  extremely  accurate,  as  the  application  of  the  method  was  the  goal,  and  not 
the  actual  results. 

In  the  fluid  component  of  the  simulations,  the  uncertain  variables  included  the  density  of  air,  the 
angle  of  attack  of  the  wing,  and  the  air  speed  of  the  plane.  For  our  model,  it  was  assumed  that 
the  aircraft  was  flying  at  an  altitude  of  9,000  meters,  which  is  safely  below  the  ceiling  cruising 
altitude  of  the  trainer.  Even  though  flight  profile  was  at  a  constant  altitude,  the  density  of  air 
could  still  have  some  variability.  At  an  altitude  of  9,000  meters,  the  nominal  density  of  air  is 
0.4671  kilograms  per  cubic  meter  [27].  A  five  percent  deviation  was  used  with  a  normal 
distribution  to  form  the  probability  distribution  for  the  air  density  input.  This  distribution  was 
chosen  because  at  Mach  0.3  there  is  a  five  percent  variation  in  the  ratio  of  the  density  of  the  fluid 
in  motion  to  the  density  of  the  fluid  at  rest  [27]. 
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Transient  flight  conditions  result  in  a  variety  of  dynamic  loads  and  conditions.  In  order  to  use  the 
vortex  lattice  method  as  a  surrogate  fluid  dynamics  model  as  described  below,  the  angle  of  attack 
needed  to  be  small  and  not  near  the  stall  angle  of  attack  since  the  method  uses  the  small  angle 
approximation  and  does  not  account  for  the  onset  of  a  stall.  Taking  these  factors  into  account, 
we  simulated  the  plane  at  a  nominal  angle  of  attack  of  eight  degrees  with  a  variability  of  plus  or 
minus  two  degrees  with  a  uniform  distribution. 

In  order  to  account  for  variation  of  bulk  air  movement  during  flight,  the  air  speed  was  also 
considered  to  be  an  uncertain  variable.  Air  speed  is  the  relative  speed  of  the  air  flowing  by  the 
plane.  As  the  wind  changes  mid  flight,  so  does  the  relative  air  speed.  In  order  for  the  vortex 
lattice  method  to  be  valid,  the  flow  field  must  be  incompressible.  For  incompressibility,  the  air 
speed  should  not  be  over  Mach  0.3.  Thus,  the  nominal  air  speed  was  set  at  91.2  meters  per 
second,  which  corresponds  to  Mach  0.3  at  an  altitude  of  9000  meters.  A  uniform  distribution 
with  a  plus  and  minus  five  meters  per  second  was  used  for  the  variation  in  airspeed. 

The  structural  surrogate  model  had  its  own  set  of  uncertain  variables.  Two  of  these  variables. 
Young’s  modulus  and  Poisson’s  ratio,  are  material  properties  of  the  7075-T6  aluminum  of  which 
most  of  the  wing  is  constructed.  The  Young’s  modulus  for  7075-T6  is  10,400  ksi,  and  Poisson’s 
ratio  is  0.33.  Not  every  sheet  of  aluminum  has  the  exact  same  mechanical  properties;  there  are 
slight  differences  in  the  heat  treatment  or  the  alloy  composition  of  two  different  sheets  of  metal. 
In  order  to  capture  this  variability,  both  the  Young’s  modulus  and  Poisson’s  ratio  were  assigned 
a  variability  of  plus  or  minus  five  percent. 

The  final  two  uncertain  variables  considered  in  the  structural  model  were  the  yield  strength  and 
the  ultimate  tensile  stress  for  7075-T6  aluminum.  There  is  a  range  of  acceptable  values  for  these 
properties,  with  the  lower  limit  set  by  the  lowest  permissible  yield  and  ultimate  strengths 
according  to  US  standards.  The  lowest  acceptable  value  for  the  yield  strength  of  7075-T6 
aluminum  is  62  ksi  [28].  However,  this  value  is  typically  as  high  as  73  ksi  [29].  The  7075-T6 
aluminum  must  have  an  ultimate  tensile  strength  of  at  least  71  ksi  [28],  but  it  can  be  as  much  as 
83  ksi  [29].  Table  1  summarizes  the  uncertain  parameters  for  the  simulations.  In  the  next  section 
the  technique  used  for  sampling  from  these  variable  distributions  is  discussed. 


Table  1.  Data  on  Fluid  and  Structure  1 

Parameters 

PARAMETER 

SYMBOL 

NOMINAL 

UNITS 

DISTRIBUTION 

VARIABILITY 

Air  density 

Pair 

0.4671 

kg/m^ 

Normai 

5% 

Angie  of  attack 

a 

8 

degrees 

Uniform 

25% 

Air  speed 

- 

91.2 

m/s 

Uniform 

5.48% 

Young's 

moduius 

E 

10400 

ksi 

Normai 

5% 

Poisson's  ratio 

v 

0.33 

- 

Normai 

5% 

Yieid  strength 

a. 

67.5 

ksi 

Normai 

7.41% 

Uitimate  strength 

O’uts 

77 

ksi 

Normai 

7.79% 
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2.23.5  Latin  Hypercube  Sampling 


The  probability  distributions  used  to  represent  the  uncertainties  in  the  fluids  and  structural 
components  need  to  be  propagated  through  the  computational  models  in  order  to  obtain  the 
variation  of  the  dependent  variables  or  system  response  quantities  (SRQs)  of  the  models.  Since 
the  uncertainties  in  the  independent  variables  were  modeled  using  probability  distributions, 
probability  theory  was  used  to  propagate  the  uncertainties  to  the  dependent  variables. 

Often,  Monte  Carlo  sampling  is  used,  in  which  random  samples  are  taken  from  the  probability 
distribution  of  each  input.  Nominal  Monte  Carlo  requires  many  samples  to  be  taken  in  order  to 
adequately  characterize  the  input  probability  distributions.  One  computer  run  must  be  done  for 
each  random  sample.  In  order  to  characterize  the  entire  set  of  distributions,  thousands  of  samples 
need  to  be  formed,  which  in  turn  requires  thousands  of  computer  runs.  This  process  is 
inefficient,  especially  if  the  tails  of  the  distribution  are  the  real  interest,  and  for  long-running 
simulation  codes,  is  actually  impossible  in  many  cases. 

A  form  of  Monte  Carlo  sampling,  known  has  Latin  Hypercube  Sampling  (LHS),  is  a  popular 
choice  of  experimental  design  when  computer  simulation  is  used  to  study  a  physical  process. 
LHS  designs  guarantee  uniform  samples  for  the  marginal  distribution  of  each  single  input.  When 
sampling  a  function  of  N  variables,  the  range  of  each  variable  is  divided  into  M  equally  probable 
intervals.  M  sample  points  are  then  placed  to  satisfy  the  Latin  hypercube  requirements;  note  that 
this  forces  the  number  of  divisions,  M,  to  be  equal  for  each  variable.  Also  note  that  this  sampling 
scheme  does  not  require  more  samples  for  more  dimensions  (variables);  this  independence  is  one 
of  the  main  advantages  of  this  sampling  scheme.  Another  advantage  is  that  random  samples  can 
be  taken  one  at  a  time,  remembering  which  samples  were  taken  so  far. 

The  LHS  process  was  applied  to  the  seven  input  parameters  and  their  probability  distributions  as 
described  in  Table  1.  Since  we  want  an  accurate  characterization  of  the  outputs,  500  samples 
were  used,  although  this  number  was  arbitrary  and  could  have  been  larger  or  smaller.  This 
means  that  for  each  input  parameter  there  were  five  samples  for  each  percentile  of  the 
distribution,  which  ensured  that  the  tails  of  the  distributions  were  well  characterized  to 
approximately  the  0.002  probability  level.  Each  sample  set  has  one  sample  member  from  each 
parameter.  The  500  sets  formed  a  matrix,  which  will  be  referred  to  as  the  LHS  sample  matrix. 
Using  a  conventional  approach,  500  simulations  of  Rocflo  and  Abaqus  would  have  been  needed, 
which  is  impractical.  However,  reducing  the  number  of  sample  sets  to  a  number  that  creates  a 
plausible  number  of  simulation  runs  on  Rocstar  (e.g.,  five  to  ten)  would  severely  limit  the 
accuracy  of  the  final  output  distributions.  The  solution  proposed  to  this  problem  of 
computational  limits  is  to  use  a  low-order  surrogate  physics  model  coupled  with  a  clustering 
methodology. 

2.2.3.6  Clustering  Methodology  Overview 

The  clustering  methodology  was  proposed  at  a  Sandia  National  Laboratories  Validation 
Workshop  [30],  and  was  initially  applied  for  uncertainty  analysis  on  CFD  in  solid  propellant 
rockets  [31].  Large  sample  sets  were  used  to  maintain  the  details  of  the  input  distributions.  Only 
a  few  high-fidelity  multi-physics  simulations  are  assumed  possible  due  to  resource  constraints. 
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Also,  it  is  assumed  that  some  of  the  LHS  sample  sets  will  produce  results  similar  to  one  another, 
and  that  one  sample  from  a  group  of  samples  that  produce  similar  results  can  be  selected  to 
represent  the  aggregate  group.  In  order  to  identify  the  representative  sample  set,  low-order 
surrogate  models  are  used  to  predict  the  approximate  response  of  the  SRQs.  These  predictions 
are  then  organized  into  clusters,  in  which  each  cluster  groups  the  predictions  that  have  similar 
results.  A  representative  sample  is  chosen  from  each  cluster  and  is  used  for  large  scale  multi¬ 
physics  simulations.  The  results  from  the  full  scale  simulations  are  interpolated  with  the  results 
from  the  surrogate  models  to  create  the  output  distribution.  This  is  a  brief  overview  of  the 
methodology  and  each  aspect  of  the  methodology  is  examined  in  more  detail  as  it  applies  to  the 
demonstration  in  the  following  sections. 

2.23.1  Surrogate  Models 

One  of  the  important  aspects  of  the  methodology  is  the  selection  of  the  low-order  surrogate 
models.  Surrogate  models  are  analytical  equations  or  fast  running  computer  simulations  that 
quickly  predict  the  SRQs.  The  surrogate  models  are  used  to  predict  the  response  of  a  SRQ  for 
each  element  in  the  LHS  matrix.  Then  the  results  produced  by  the  surrogate  models  are 
clustered  as  describe  in  section  2. 2.3. 6,  and  a  representative  sample  set  of  input  values  is  chosen 
for  each  cluster.  The  main  purpose  of  the  surrogate  model  is  to  define  the  shapes  of  the  SRQ 
distributions  in  relation  to  the  input  probability  distributions.  The  surrogate  models  capture  the 
trends,  not  the  exact  values,  of  the  output  parameter  correctly.  The  models  chosen  should  give 
similar  result  trends  to  those  that  would  be  seen  if  the  Rocstar/Abaqus  models  were  run  for  each 
LHS  sample  member.  A  few  full-scale  detail  simulations  are  then  performed  to  allow 
interpolating  the  surrogate  results  into  a  more  “accurate”  space. 

Since  these  simulations  contain  a  fluid  and  structure  component,  surrogate  models  were 
developed  that  captured  the  uncertainty  trends  in  both.  The  main  focus  of  the  project  was 
predicting  the  structural  response  rather  than  the  fluid  response;  consequently,  several  structural 
surrogate  models  of  differing  degrees  of  fidelity  were  tested  to  investigate  the  impact  of  the 
surrogate  models  on  the  output  parameter’s  probability  distribution. 

2.2.3.7.1  Fluids  Surrogate  Model 

For  the  fluids  component  of  the  problem,  the  surrogate  model  involved  the  three  uncertain 
variables:  the  angle  of  attack,  air  density,  and  the  airspeed.  These  variables  were  propagated  into 
the  model  in  the  form  of  lift  on  the  wing.  The  total  lift  on  an  airplane  can  be  found  using 
Equation  1: 


(16) 

where  Cl  is  the  lift  coefficient,  pair  is  the  air  density,  and  Uoo  is  the  freestream  air  velocity.  The 
angle  of  attack  is  incorporated  into  the  lift  coefficient,  and  thus  is  also  a  variable  of  the  total  lift. 
S  represents  the  planform  area  of  the  wing;  however  this  was  not  an  uncertain  variable  since  the 
geometry  of  the  plane  was  assumed  constant  within  tolerances  that  will  not  affect  the 
aerodynamics  of  the  aircraft. 
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The  lift  coefficient  and  thus  the  total  lift  were  calculated  using  the  Matlab  program  Tornado  [32]. 
This  program  uses  the  vortex  lattice  method  to  describe  the  airflow  around  the  wing  of  an 
aircraft.  A  few  assumptions  need  to  be  made  about  the  flow  field  in  order  for  this  method  to  be 
valid.  The  flow  field  must  be  incompressible  and  irrotational.  Also,  viscous  forces  are  ignored 
and  the  angle  of  attack  must  be  small,  so  that  the  small  angle  approximation  may  be  used. 

In  Tornado,  the  geometry  of  the  wing  of  the  trainer  had  to  be  approximated.  Inputs  included  the 
root  chord  length,  the  dihedral  and  sweep  of  the  wing,  the  taper  ratio,  and  span  of  the  wing.  In 
order  to  create  the  correct  airfoil,  coordinates  of  the  NACA  65A004.8  airfoil  were  normalized 
and  then  imported  into  Tornado.  The  geometry  and  airfoil  of  the  trainer  represent  constant 
properties  that  do  not  contain  uncertain  variables.  The  uncertain  variables  appeared  in  the 
description  of  the  state,  or  environment,  of  the  aircraft. 

The  fluid  uncertainty  variables  were  propagated  through  Tornado  and  the  total  lift  was  recorded. 
The  total  lift  was  used  in  the  structural  surrogate  models  to  determine  the  magnitude  of  the  lift 
distribution. 


2.23.1.2  Structural  Surrogate  Models 

Each  of  the  structural  surrogate  models  was  created  using  the  same  basic  assumptions.  The  wing 
was  considered  to  be  entirely  made  from  A1  7075-T6.  Next,  the  wing  was  assumed  to  be  under 
an  elliptical  lift  distribution,  as  described  in  Equation  2,  where  E  is  the  length  of  the  wing  and  x 
is  the  distance  along  the  wing  measured  from  the  wing  root.  The  variable  /o  refers  to  the 
magnitude  and  is  calculated  in  Equation  3. 


( 

F(x)  =  f,  1 

V 


/o 


(Lift) 


(17) 

(18) 


It  was  assumed  that  the  wing  was  under  pure  bending,  i.e.,  the  wing  to  behaved  like  a  beam.  The 
wing  was  approximated  as  a  cantilever  beam  since  the  dimensions  of  its  airfoil  were  small 
compared  to  the  longitudinal  axis  of  the  wing.  Although  the  cross  section  of  the  wing  does  not 
change  in  shape,  there  is  a  gradual  decrease  in  size.  Because  of  this,  the  moment  of  inertia  was 
averaged  in  order  to  make  it  constant  and  simplify  the  calculations.  This  average  inertia  was 
taken  at  the  Mean  Aerodynamic  Chord  (MAC).  The  MAC  is  the  characteristic  chord  at  which  the 
average  properties  of  the  airfoil  can  be  determined.  When  the  airfoil  at  this  chord  was  modeled 
in  Abaqus,  a  moment  of  inertia  of  16719.3  in"^  was  calculated.  This  value  was  assumed  to  be  the 
average  inertia  of  the  airfoil  across  the  entire  wing.  The  cross  section  contained  an  effective  skin 
thickness  that  incorporated  the  top  and  bottom  flanges  of  seven  main  spars.  Drag  was  ignored 
and  the  twisting  of  the  wing  due  to  the  shear  stress  was  not  modeled. 
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2.23.1. ?>  Castigliano’s  Theorem  Model 


The  first  structural  surrogate  model  uses  Castigliano’s  theorem,  i.e.,  the  derivative  of  the  total 
strain  energy  of  an  elastic  system  with  respect  to  the  force  acting  at  a  specific  location  is  equal  to 
the  displacement  at  that  location. 


1  rM" 

U  =  r  ^  dx 

(19) 

2{eI 

dU 

(20) 

dQi 

Equation  4  shows  how  the  strain  energy  (U)  for  the  idealized  wing  is  calculated.  In  order  to 
simplify  this  analysis,  it  was  assumed  that  Young’s  modulus,  E,  and  the  moment  of  inertia,  I,  of 
the  airfoil  were  both  constants  for  any  single  realization  of  the  model.  To  find  the  moment 
associated  with  the  wing,  Mz,  it  was  assumed  that  an  elliptic  lift  curve  was  distributed  over  the 
entire  span  of  the  wing.  In  Equation  5  the  term  qi  refers  to  the  displacement  at  a  specific  point  of 
the  beam,  while  Qi  is  the  generalized  force  at  that  point.  Since  we  were  interested  in  the 
displacement  at  the  wing  tip,  a  “dummy  load”  (Pd)  was  required  at  the  wing  tip.  In  order  to 
account  for  this  dummy  load.  Equation  5  was  adjusted: 


w{L)  = 


lim  (du'^ 
1  D  ^  Z)  J 


(21) 


where  w(E)  is  the  displacement  of  the  wing  tip  in  the  z  direction.  Each  sample  member  in  the 
EHS  sample  was  processed  with  this  model  (using  the  fluids  surrogate  model  from  section 
2. 2. 3.7.1  to  produce  the  load),  and  an  estimate  of  the  wing  tip  displacement  was  generated  for 
each  sample. 

2.2.3.7.4  Euler-Bernoulli  Beam  Model 


Euler-Bemoulli  beam  theory  was  the  second  surrogate  model.  Equation  7  shows  the  wingtip 
displacement,  where  /o  is  the  load  factor,  E  is  the  length  of  the  wing  or  half  the  span,  E  is 
Young’s  modulus,  and  I  is  the  second  area  moment  of  inertia. 


19/oL^ 

360E/ 


(22) 


The  load  factor, /o,  was  determined  by  the  total  lift  which  was  calculated  using  the  fluid  surrogate 
model.  The  uncertainty  is  propagated  through  Young’s  modulus  and  the  load  factor.  The 
surrogate  model  was  run  500  times  with  the  EHS  matrix  and  the  displacement  recorded. 
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2.2.3.7.5  Abaqus  Beam  Model 


The  Abaqus  wing-beam  model  used  the  Meshed  Beam  Cross-Section  (MBCS)  approach 
described  in  the  Abaqus  documentation  [33].  MBCS  allows  structures  with  complex  geometry, 
which  exhibit  characteristics  of  beams,  to  be  modeled  without  having  to  create  a  3D  model.  It  is 
a  two-step  process  in  which  the  cross  sectional  properties  are  calculated  separately,  in  meshed 
cross  section  models,  and  then  assigned  to  the  corresponding  segments  of  the  beam  model  for  the 
beam  analysis.  This  surrogate  model  had  the  highest  degree  of  fidelity  of  the  three  solid 
mechanics  surrogates  because  it  accounted  for  varying  geometry. 

In  order  to  model  the  swept  and  tapered  wing,  the  wing  was  partitioned  into  five  segments  that 
were  separated  by  ribs.  Each  segment  had  a  constant  cross  section  that  was  defined  by  the  cross 
section  of  the  MAC  for  that  segment.  This  resulted  in  decreasing  cross  sections  that  modeled  the 
geometry  of  the  wing.  Each  cross  section  modeled  the  skin  and  spars  as  one  piece.  MBCS  does 
not  model  the  contact  or  riveting  between  structural  components.  The  skin-spar  part  was 
modeled  using  warp  elements  and  assumed  linear  isotropic  elastic  material.  This  material 
captured  the  uncertain  parameters  of  Young’s  modulus  and  Poisson’s  ratio.  From  the  cross 
section  models,  an  out-of-plane  warping  function  and  the  cross  section  properties  were 
calculated  (area  moment  of  inertia,  torsion  constant,  shear  center,  etc.).  These  were  utilized  in  the 
beam  analysis. 

Using  MBCS,  the  wing  was  modeled  as  a  cantilever  beam  with  linear  beam  elements.  These 
elements  used  Timoshenko  beam  theory.  The  cross  sectional  properties  that  were  defined  in  the 
meshed  cross  section  models  were  assigned  to  the  linear  beam  elements  that  represent  that 
segment  of  the  wing.  In  Abaqus,  a  pressure  distribution  cannot  be  applied  to  linear  beam 
elements;  consequently,  for  each  segment  a  force  resultant  was  applied.  The  beam  analysis  was 
conducted  and  the  wing  tip  displacement  found  for  each  of  the  EHS  sample  members. 

2.2.3.8  Clustering 

Restrictions  on  the  number  of  computer  simulations  that  can  be  performed  occur  frequently.  In 
our  case,  the  time  needed  to  run  the  fluid  flow  solver,  Rocflo,  was  the  main  restriction  to  the 
number  of  simulations  that  could  be  performed.  We  did  not  have  time  to  run  500  separate  Rocflo 
simulations,  so  we  limited  the  number  of  full  scale  runs  to  six.  To  select  the  six  EHS  sets,  the 
EHS  matrix  was  processed  with  the  fluid  surrogate  model.  The  resulting  wing  loads  were  then 
used  in  each  structural  surrogate  model  and  the  wing  tip  displacements  found  for  each  sample 
member  in  the  EHS  matrix.  In  this  methodology,  the  number  of  clusters  required  is  one  less  than 
the  total  number  of  high-fidelity  Rocstar/Abaqus  simulations  that  are  planned. 

Using  the  wing  tip  displacements.  The  500  displacements  for  each  surrogate  model  were 
clustered  into  five  groups.  Each  cluster  represents  a  set  of  samples  that  predicted  similar 
displacements  of  the  wing  tip.  This  clustering  process  was  done  using  SPSS,  a  statistical 
package  that  employs  a  K-means  clustering  algorithm.  This  type  of  algorithm  is  a  partitional 
method  where  the  number  of  clusters  must  be  predefined  by  the  user.  The  clusters  for  the 
Castigliano’s  Theorem,  the  Euler-Bemoulli  Beam  Theorem,  and  tht  Abaqus  Beam  surrogates  are 
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shown  in  Figures  7  to  9.  Each  circle  on  the  graph  represents  the  wing  tip  displacement  of  one  of 
the  members  of  the  LHS  matrix. 

The  representative  samples  that  were  to  be  run  in  the  full  scale  multi-physics  simulations  were 
selected  based  on  these  clusters.  There  are  different  ways  in  which  this  could  be  done  [34].  The 
effects  of  the  extreme  values  are  of  importance  since  the  goal  was  to  create  probability 
distributions  for  various  properties  of  the  wing.  The  best  way  to  capture  the  extreme  values  is  to 
take  the  first  and  last  sample  member  of  each  cluster  after  they  are  sorted  in  ascending  order  by 
the  wing  tip  distribution.  Using  this  method,  a  total  of  ten  sets  of  input  parameters  were 
collected  from  the  five  clusters. 


However,  the  lowest  sample  of  a  cluster  and  the  highest  sample  of  a  previous  cluster  have  almost 
identical  responses  since  the  range  of  wing  tip  displacements  was  almost  continuous,  as  seen  in 
Figures  6  to  8.  Since  the  Rocstar  responses  should  be  similar,  it  is  only  necessary  to  take  one  of 
the  two  values  from  each  cluster  interface.  Thus,  by  taking  only  the  sample  member  with  the 
highest  displacement  for  each  cluster  and  also  the  lowest  sample  in  the  lowest  cluster,  the 
number  of  simulations  was  reduced  to  six. 

The  CFD  parameters  from  the  six  selected  sample  members  were  used  in  six  Rocflo  runs.  Each 
of  those  Rocflo  runs  was  performed  with  its  unique  set  of  parameters,  and  the  pressures  that 
developed  on  the  surfaces  of  the  wing  during  the  simulation  were  extracted  and  then  applied  to 
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an  Abaqus  model  using  the  corresponding  Solid  Mechanics  parameters  for  the  given  sample.  The 
Abaqus  model  was  then  run,  and  a  wing  tip  displacement  was  found.  This  provides  “high- 
fidelity”  estimates  of  the  wing-tip  displacement  for  six  cases.  These  six  high-fidelity  results  were 
then  mapped  onto  the  distribution  estimated  from  the  surrogate  model  results.  This  is  achieved 
by  interpolating  between  the  high-fidelity  results  based  on  the  surrogate-model  results.  This 
effectively  “pinned”  the  SRQ  distribution  at  six  points  (including  the  extrema)  based  on  the  high- 
fidelity  results,  and  then  the  shape  of  the  distribution  between  those  points  was  based  on  the  low- 
order  surrogate  results.  Thus,  the  accuracy  of  the  surrogate  models  is  not  important;  they  just 
need  to  predict  the  trends  correctly. 


Figure  8:  Cluster  ranges  for  Euler-Bernoulli  model 
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Figure  9:  Cluster  ranges  for  the  Abaqus  beam  model 


2.2.3.9  Results 

One  of  the  goals  of  this  projeet  was  to  investigate  the  effects  that  different  surrogate  models 
would  have  on  the  final  output  probability  distributions.  The  trends  of  the  output  distributions 
are  most  impacted  by  the  surrogate  models.  In  Figure  9,  the  normalized  wing  tip  displacements 
for  each  surrogate  model  are  shown  for  each  sample  set  in  the  LHS  matrix.  The  results  trends 
are  almost  exactly  the  same.  This  is  to  be  expected  because  Castigliano’s  Theorem  and  the 
Euler-Bemoulli  surrogate  models  are  different  methods  to  reach  the  same  conclusion.  The 
Abaqus  beam  model  provided  more  accurate  values  for  the  displacement;  however,  because  the 
analysis  was  linear,  it  provided  a  very  similar  trend  to  the  other  models. 
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Figure  10:  The  normalized  displacements  comparing  trends  of  the  surrogate  models 


The  final  step  in  the  uncertainty  assembly  has  been  attempted  using  this  methodology.  The 
results,  however,  were  disappointing  from  the  perspective  that  the  displacements  calculated  in 
Abaqus  based  on  the  transferred  CFD  pressures  were  very  small.  This  may  be  due  to  limited 
runtime  for  the  CFD  simulations,  the  limited  angle  of  attack  able  to  be  simulated  due  to  the 
surrogate  model  limitations,  or  errors  in  the  pressure  transfer  algorithms  that  were  developed  for 
this  work.  The  direction  of  this  project  was  redirected  before  it  was  possible  to  ascertain  where 
the  problem  exists.  Thus,  while  the  overall  process  appears  to  be  feasible  in  a  Fluid-Structure- 
interaction  problem,  a  complete  end-to-end  assessment  will  be  a  future  project. 
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3.0  CONCLUSIONS  AND  RECOMMENDATIONS 


Probabilistic  methods  applied  to  ensuring  structural  integrity  show  promise  as  they  have  for  the 
past  30  years.  The  challenge  does  not  seem  to  be  in  developing  new  more  sophisticated 
probabilistic  techniques,  but  in  getting  the  techniques  into  common  practice.  It  needs  to  become 
relatively  easy  for  the  average  engineer  without  a  lot  of  expertise  to  use  probabilistic  methods. 
The  way  to  do  this  is  to  make  incremental  changes  to  existing  design  and  maintenance  processes 
by  enabling  engineers  to  readily  obtain  more  and  better  information  upon  which  to  base  any 
decisions  about  the  airworthiness  of  an  aircraft  upon.  This  should  be  the  focus  of  future  efforts 
on  Risk-Quantified  Structural  Design  and  Evaluation. 
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